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Using the nonequilibrium Green function formalism combined with density functional theory, 
we study finite-bias quantum transport in Ni|Gr n |Ni vertical heterostructures where n graphene 
layers are sandwiched between two semi-infinite Ni(lll) electrodes. We find that recently predicted 
"pessimistic" magnetoresistance of 100% for n > 5 junctions at zero bias voltage Vb — ► 0, persists 
up to Vb ~ 0.4 V, which makes such devices promising for spin-torque-based device applications. In 
addition, for parallel orientations of the Ni magnetizations, the n = 5 junction exhibits a pronounced 
negative differential resistance as the bias voltage is increased from V& = V to Vb — 0.5 V. We 
confirm that both of these nonequilibrium transport effects hold for different types of bonding of 
Gr on the Ni(lll) surface while maintaining Bernal stacking between individual Gr layers. 

PACS numbers: 72.25.Mk, 73.43.Qt, 75.47.-m, 72.80.Vp 



I. INTRODUCTION 

A magnetic tunnel junction (MTJ) consists of an ultra- 
thin insulating barrier which separates two metallic fer- 
romagnetic (F) layers with variable magnetization direc- 
tion. The MTJs based on transition metals or their alloys 
and an epitaxial MgO barrier 1 are the present workhorse 
of both commercial and basic research spintronics. For 
example, MgO-based MTJs are the core elements of read 
heads in hard drives or in magnetic random-access mem- 
ory devices 2 that are operated by the current-induced 
spin-transfer torque (STT). In the STT phenomenon, 
spin current of large enough density injected into a free 
F layer either switches its magnetization from one static 
configuration to another or generates a dynamical situa- 
tion with steady-state precessing magnetization. 3 Thus, 
the usage of MTJs in STT-based spintronic devices ne- 
cessitates 2 ^ a compromise between large current density 
(which requires low junction resistance to avoid damage) 
driven by finite bias voltage and readability (which re- 
quires large magnetoresistance). 

One of the great successes of first-principles quan- 
tum transport modeling has been a conjecture^ of very 
large "optimistic" tunneling magnetoresistance, TMR = 
(G P - G A p)/G A p x 100% > 1000%, in Fe|MgO|Fe(100) 
MTJs where Gp (Gap) is conductance for parallel (an- 
tiparallel) orientation of the Fe magnetizations. This 
prediction has ignited intense fabrication efforts reach- 
ing TMR of about 200 % at room temperature 5 which, 
although undoubtedly correlated with the crystallinity of 
MgO barrier, is difficult to reconcile with first-principles 
predictions. 1 The origin of the discrepancy is the sensi- 
tivity of spin injection and TMR to details of difficult- 
to-control interfacial disorder as revealed by a number 
of theoretically investigated scenarios (such as the in- 
termixing of Fe and MgO, 6 oxygen vacancies at or near 
the Fe|MgO interface, 7 or substoichiometric FeO layers 



with small oxygen concentrations^). In addition, TMR 
in MgO-based MTJs decays precipitously 519110 ! with in- 
creasing bias voltage where the specific features of the 
decay are also sensitive to the type of interfacial disor- 
der P 

These issues could be resolved by searching for new 
material systems which would ensure perfect spin fil- 
tering in the absence of disorder while being much less 
sensitive to the presence of interfacial disorder in real- 
istic junctions. For example, the recent first-principles 
analysi d 11 * 12 ! has brought an example of such system — 
Ni|Gr n |Ni junctions — where n layers of graphene (Gr n ) 
are sandwiched between two Ni electrodes as illus- 
trated in Fig. [TJ Graphene is recently discovered^ two- 
dimensional (2D) allotrope of carbon where electronic 
states of a single layer Gri or multilayers Gr n close to 
the Fermi energy are located around the high symmetry 
K point in reciprocal space. The Ni|Gr n |Ni junction ex- 
ploits very small mismatch of 1.3% between the in-plane 
lattice constant of Gr and the surface lattice constant 
of Ni(lll), as well as the fact that majority spin states 
of Ni are absent in a large region around the K point. 
These two features combined lead to perfect spin filtering 
for n > 5, as quantified by the "pessimistic" magnetore- 
sistance MR = (G P - G A p)/G P x 100% « 100 % (the 
"optimistic" MR diverges since Gap vanishes for large 
n). The three times smaller lateral lattice mismatch com- 
pared to the 3.8% for conventional Fe|MgO|Fe junctions 
should also reduce some of the strain and amount of de- 
fects that otherwise limit the thickness and degrade the 
efficiency of spin injection in MgO-based MTJs. 

However, very little is known about nonequilibrium 
transport driven by finite bias voltage Vb in Ni|Gr n |Ni 
junctions. This is partly due to the fact that stan- 
dard first-principles electronic transport tools employed 
to capture electronic and magnetic structure at inter- 
faces, such as layer Korringa-Kohn-Rostoker approach 
applied 1 to MgO-based MTJs or tight-binding muffin 
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FIG. 1. (a) Schematic view of Ni|Grs|Ni junction where Grs 
represents five layers of graphene and Ni is (111) fee nickel. 
The device extends to infinity along the transverse directions 
while Ni electrode (orange) are semi-infinite in the longitu- 
dinal (transport) direction. The two investigated types of 



bonding 



'for Gr on the Ni(lll) surface are illustrated 



in panel (b) , as AB configuration where the two carbon atoms 
in the graphene unit cell cover Ni atoms in layers A (surface) 
and B (second layer), and panel (c) as AC configuration in 
which carbon atoms are placed directly above the Ni atoms 
in layers A (surface) and C (third layer). Here ABC refers to 
three close-packed layers within a fee crystal. 



tin orbital wave-function matching scheme appliecP^EU 
to Ni|Gr n |Ni junctions, become very cumbersome^ to 
use at finite Vb where one has to compute the charge 
redistribution^ due to current flow by evaluating the 
nonequilibrium density matrix p — a procedure which en- 
sures the gauge invariance^ of the current- volt age I- V 
characteristics. The nonequilibrium Green function for- 
malis m combi ned with density functional theory (NEGF- 
DFT) j 15 ' 17 ' 1 ^ where DFT part of the calculations is im- 
plemented in the basis of local orbit als, makes it rela- 
tively straightforward to obtain p. 

Here we show how to use efficiently spin- and kn- 
resolved NEGF-DFT framework to understand nonequi- 
librium transport through to Ni|Gr n |Ni junctions de- 
picted in Fig. [I] for parallel (P) or antiparallel (AP) ori- 
entation of the Ni magnetizations. Our principal results 
are shown in Fig. [2] and Fig. [3] In Fig. [2|a), we first 
confirm the result of Ref. [11] and [12] about the zero bias 
"pessimistic" MR reaching 100% for barriers composed of 
n > 5 graphene layers and, moreover, in Fig.[2|b) we pre- 
dict that such maximized MR would persist even at finite 
Vb i$ 0-4 V. Figure [2^b) also suggests that bias voltage 
dependence of MR can be employed experimentally to 
determine the type of bonding configuration [illustrated 
in Figs, [ljb) and (c)] for Gr on the Ni(lll) surface. 

Furthermore, Fig. [3jc) shows that Ni|Gr n |Ni junction 
with P orientation of the Ni magnetizations will exhibit 
pronounced negative differential resistance (NDR), where 
total charge current first increases and then decreases as 
the bias voltage is increased from V& = V to Vb — 0.5 V 
(or symmetrically in the opposite direction). The origin 
of NDR is explained in Fig. [5] by plotting the position- 
dependent local density of states (LDOS) across the junc- 
tion. 
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FIG. 2. (a) The "pessimistic" TMR for Ni|Gr n |Ni junctions 
as a function of the number of graphene layers n and for two 
different, AB and AC, bonding configurations for Gr on the 
Ni(lll) surface illustrated in Figs.[TJb) and[TJc), respectively, 
(b) The "pessimistic" TMR for n — 5 junction versus finite 
bias voltage for AC and AB bonding configurations. 



The paper is organized as follows. In Sec.[ll]we discuss 
vertical Ni|Gr n |Ni heterostructure setup and how to tune 
the spin- and ky -resolved NEGF-DFT framework in or- 
der to describe properly magnetism around its interfaces. 



Section III discusses magnetoresistance at finite bias volt- 
age, as well as the unusual conduction properties of gen- 
eral vertical graphene heterostructures, whose fabrica- 
tion has been initiated recently^ which make Ni|Gr n |Ni 
junctions different from either conventional MTJs or spin 
valves. In Sec. Ill, we discuss NDR in Ni|Gr n |Ni junc- 



tions driven by finite bias voltage. We conclude in Sec.[V| 



II. THE VERTICAL HETEROSTRUCTURE 
SETUP AND TUNING OF NEGF-DFT 
FRAMEWORK FOR ITS MODELING 

The disorder- free junction shown in Fig. [IJa) consists 
of up to seven graphene layers arranged in Bernal stack- 
ing^ which serve as the barrier separating the two semi- 
infinite Ni electrodes. The junction is infinite in the 
transverse direction, so that its transverse periodicity 
requires /c-point sampling. 9 The spin injection and spin 
filtering in ferromagnetic multilayers depends not only 
on the properties of the F electrodes but also on geome- 
try, bonding and electronic and magnetic structure of the 
contact region, as emphasized by the studies^ Q f MgO- 
based MTJs. Therefore, we consider two different Gr on 
the Ni(lll) surface bonding configurations illustrated in 
Figs. Qb) and [TJc). 

We note that DFT calculations employing different 
approximations for the exchange-correlation functional 
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FIG. 3. The finite-bias transmission function T min (E,V b ) 
for Ni|Gr n |Ni junction in AC bonding configuration at the 
Ni(lll)|Gr interface [Fig. [IJc)] for (a) P and (b) AP ori- 
entations of the Ni magnetizations. Since in P orientation 
minority spin contribution dominates, while in AP setup 
both minority and majority spins contribute the same, only 
T min (E, Vb) is presented here for both P and AP orientations 
with curves at different Vb shifted along the y-axis for clar- 
ity. Panels (c) and (d) show /- V characteristics for P and 
AP orientation, respectively. The NDR is conspicuous in P 
orientation in panel (c) for both AC and AB bonding config- 
urations. 



(such as local density approximationpP^ generalized 
gradient appro ximation^ and van der Waals density 
functiona l 21 * 22 *) have yielded contradictory conclusions 
about the AC bonding configuration being the most sta- 
ble energetically and the corresponding binding distance. 
The recent random phase approximation (RPA) calcula- 
tions^ have resolved this controversy and demonstrated 
the conflicting results are due to a delicate interplay be- 
tween covalent and dispersive interactions which is not 
captured by the DFT functionals. Also, the scanning 
tunneling microscopy imaging 2 -^ shows that perfectly or- 
dered epitaxial graphene layers can be prepared by ele- 
vated temperature decomposition of hydrocarbons where 
domains are larger than the terraces of the underlying 
Ni(lll) surface. 

T he N EGF-DFT framework was originally devel- 
opecP^HI to treat quantum transport through small 
molecules attached to metallic electrodes. Its applica- 
tion to modeling of charge and spin transport in MTJs re- 
quires careful tuning of pseudopotentials and basis sets in 
order to obtain an accurate description of the band struc- 
ture near the Fermi level which is particularly important 
for the studies of spin-polarized transport. For example, 
pseudopotentials and localized basis sets that reproduce 
the electronic structure of the F electrode and barrier 
material alone do not necessarily reproduce the electronic 
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FIG. 4. The comparison of first-principles computed band 
structure of a periodic • • • Ni|Gr5|Ni|Gr5 • • • superlattice [with 
AC bonding configuration for Gr on the Ni(lll) surface] ob- 
tained using either real-space grid PAW method implemented 
via the GPAW code 25 or SZP basis of localized orbitals on C 
atoms and DZP basis on Ni atoms, together with pseudopo- 
tentials, implemented via the ATK code" 



structure of the more complicated F-electrode| barrier in- 
terfacesP 

In order to capture accurately the electronic and 
magnetic structure around interfaces, we first com- 
pute the band structure of a periodic superlattice 
• • • Ni|Gr 5 |Ni|Gr 5 • • • using DFT based on the projector 
augmented wave (PAW) methodology and wave func- 
tion representation on uniform real-space grids as im- 
plemented in the GPAW code 25 where we choose grid 
spacing of 0.18 A. Within the PAW formalism one 
works implicitly with the all-electron wave functions and 
has access to the (frozen) core states which make the 
method applicable to a broad range of systems (includ- 
ing materials with strongly localized d or / electrons 
that can be problematic to describe with pseudopoten- 
tials). The same band structure is then recomputed us- 
ing DFT where wavefunctions are represented in terms 
of the linear combination of atomic orbitals (LCAO) and 
the behavior of the core electrons is described by norm- 
conserving Troullier-Martins pseudopotentials, as imple- 
mented in the ATK code. 26 In the ATK-based calcu- 
lations, we choose single ("-polarized (SZP) basis on C 
atoms and double ("-polarized (DZP) basis on Ni atoms. 
The Brioullin zone of the superlattice was sampled by 
12 x 12 x 100 /c-point grid, and the charge density and 
potentials were determined on a real-space grid with a 
mesh cutoff energy of 150 Ry. This was sufficient to 
achieve a total energy convergence of better than 0.01 
meV/unit cell in the self-consistent loop. 

The excellent agreement we achieve in Fig. [4] between 
the real-space grid PAW and LCAO pseudopotential 
DFT calculations, where Perdew-Burke-Ernzerhof (PBE) 
parameterization of the spin-polarized generalized gradi- 
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ent approximation (GGA) for exchange-correlation func- 
tional has been used in both cases, also selects the correct 
parameters to be used for LCAO pseudopotential part 
of NEGF-DFT analysis of two-terminal Ni|Gr n |Ni junc- 
tions discussed below. The active region of the device in 
Fig. fija) simulated by the NEGF-DFT code consists of 
7 Ni(lll) layers on the left, n layers of Gr, and 6 lay- 
ers of Ni(lll) on the right. This active region is first 
relaxed until the maximum force component goes below 
0.02 eV/A per atom, and then attached to two semi- 
infinite ideal Ni electrodes. 

The NEGF formalism for steady-state transport op- 
erates with two central quantities — the retarded G(E) 
and the lesser Green functions G < (E) — which describe 
the density of available quantum states and how elec- 
trons occupy th ose st ates, respectively. In the NEGF- 
DFT framework j 15 * 17 * 1 ^ the Hamiltonian is not known in 
advance and has to be computed by finding the converged 
spatial profile of charge via the self-consistent DFT loop 
for the density matrix p = ^ J dEG < {E) whose diag- 
onal elements give charge density. 15 The Hamiltonian 
matrix H in the local orbital basis {^} is composed of 
elements = (^|^Ks|0j) 5 where #ks is the effective 
Kohn-Sham Hamiltonian obtained from the DFT self- 
consistent loop and the overlap matrix S has elements 

In the coherent (i.e., in the absence of electron- 
phonon or electron-electron dephasing processes) trans- 
port regime, only the retarded Green function 

Gk„ =[£S-H£,,-S£, k|| -S£ ik|| ]-\ (1) 

of the active device region is required to post-process the 
result of the DFT loop by expressing the current between 
the left (L) and the right (R) electrodes 

I a (V b ) = I dkJdET«(k b E,V b )[f(E-VL)-f(E- m )}. 
Jbz J 

(2) 

The electrodes are assumed to be attached to macro- 
scopic reservoirs at infinity characterized by the Fermi 
function f(E — ^l,r)^ so that bias voltage driving the 
nonequilibrium transport is given by hl — I^R — eV b . 
Here we resolve all quantities in minority and majority 
spin channels — a = min, maj — assuming irrelevance of 
spin-orbit coupling or spin-flip scattering. The spin- and 
k|| -resolved transmission function of coherent transport 
is given by 

T CT (k,|, E, V b ) = Tr {r£, k|| (E)G^ (£7)r£, k|| (E)G^(E)} , 

(3) 

where the level broadening matrices 
r L(K) )k|| (^)=^ S L(K) )k|| (^)- S Lt), k|| (^)] are ex- 
pressed in terms of the retarded self-energies ^l(r) k y (^) 
of semi-infinite ideal Ni electrodes. In order to con- 
verge integration over the (conserved in the absence 
of disorder) transverse wavevector ky in Eq. ([2|, we 
find it necessary to use a dense grid 301 x 301 of 
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TABLE I. The approximative values for the linear-response 
conductances, in units of Q _1 (/im) -2 , for Ni|Gr5|Ni junctions 
in AC [Fig.JIJc)] and AB [Fig.[TJb)] bonding configuration for 
Gr on the Ni(lll) surface and for P and AP orientations of 
the Ni magnetizations. The third row shows the same con- 
ductances computed in Ref. 1 for Fe|MgO|Fe(100) MTJ con- 
taining six-layer MgO barrier. 



^-points in the corresponding 2D BZ. This proce- 
dure yields the bias-dependent transmission function 
T(E,V b ) = J BZ dk\\T(lL\\,E,V b ) plotted in Figs.p^a) and 
lb). 



III. MAGNETORESISTANCE AT FINITE BIAS 
VOLTAGE 

The Ni|Gr n |Ni multilayered heterostructure is not 
a conventional MTJ. Unlike MgO-based MTJs where 
linear-response (V b — >• 0) conductances Gp m = I min /V b 
and Gp aj = I ma i /V b decay exponentially^ with increas- 
ing number of MgO layers, in the case of Ni|Gr n |Ni junc- 
tion Gp m is independent of n for n > 4 (apart fro m an 
even-odd oscillation as a function of the thickness n)! 11 * 12 * 
On the other hand, Gr n acts as a tunnel barrier for ma- 
jority spin electrons causing Gp aj to decay exponentially 
with n. The spin-resolved linear-response conductances 
for Ni|Gr5 1 Ni junctions are compared in [I] with the same 
conductances^ for Fe|MgO|Fe MTJ containing MgO bar- 
rier of similar thickness as our Gr$ barrier. 

The recent first-principles analysis^ of different 
metal | Gr n | metal junctions for n < 4, assuming reason- 
able metal-graphene epitaxial relationships, has delin- 
eated conditions for Gr n to behave effectively as a tun- 
nel barrier causing exponential decay of the conductance 
with increasing n which requires crystal momentum mis- 
match between the bulk Fermi-level states in the metal- 
lic electrode and those in the Gr n barrier. Furthermore, 
the recent experiments^ measuring /- V characteristics 
of metal | Gr n | metal vertical junctions (with Ti/Pt used 
as top and bottom metal electrodes) has demonstrated 
application of the bias voltage up to \V b \ < 1 V with- 
out encountering catastrophic breakdown while showing 
transitions from Ohmic I (xV b (at very low bias) to power 
law / ex V™ (m > 2) characteristics. 

In conventional MTJs, tunneling rates are higher if 
there are similar or identical states on both sides of the 
barrier. Therefore, the tunneling electrons need not only 
to get through the barrier but there must be a state of the 
correct symmetry on the other side to accept themP^ This 
effect is part of the reason for the commonly observed 
decrease in TMR with V b since as the bias increases the 
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FIG. 5. The position-dependent LDOS from left to right electrodes in Ni|Gr5|Ni junction, in AC bonding configuration at the 
Ni(lll)|Gr interface and P orientation of the Ni magnetizations, at different bias voltages V b . The electrochemical potentials 
jj/L and jiR of the two Ni electrodes are marked by dashed horizontal lines while the zero of energy is set at (/jll + hr)/2. The 
LDOS exhibits high values in the Ni electrodes (white regions), while the central colored region corresponds to the Grs barrier. 
The dashed ovals indicate the position of the resonant states which contribute to transport. Note that a strong coupling of the 
resonant states of the electrodes and the Gr barrier at a given energy level is required for large transmission T(E, Vb) through 
the junction. 



states on o pposite sides of the barrier for P orientation 
differ moreP^™ On the other hand, "pessimistic" MR 
in Ni|Gr n |Ni remains 100% up to V b < 0.4 V for Gr bar- 
rier of thickness n > 5, as shown in Fig. [2|b). 



IV. NEGATIVE DIFFERENTIAL RESISTANCE 

Figures |3jc) and (3^d) plot the I- V characteristics for 
P and AP orientations of the Ni magnetizations where 
the total charge current is / = 7 mm + I m& i . Since in the 
AP orientation the bias-dependent transmission T(E, Vb) 
in Fig. |3jb) is nearly flat around the Fermi level, the 
/- V characteristics in Fig. [3jd) is linear up to the volt- 
age Vb ~ ±0.6 V. However, in P orientation the total 
charge current / sharply increases to reach its maxi- 
mum value at Vb ~ ±0.12 V and then drops, thereby 
exhibiting a pronounced NDR. This feature can be ex- 
plained using T(E,Vb) curves plotted for AC configura- 
tion in Fig.[3ja). At lower the transmission resonance 
(around E — Ep — 0.0 eV) falling into the bias window 
(marked by dashed wedge) contributes to the peak in the 
I-V characteristics. However, this resonance gets dimin- 
ished with increasing Vb which eventually shuts off the 
current flow when Vb ~ 0.5 V is reached. The current 
is allowed to flow again when the new resonance around 
E — Ep « -0.45 eV enters the bias window V b « -0.7 V. 

Further insight into the microscopic mechanism behind 
NDR in P orientation of magnetizations in Ni|Gr n |Ni 
junctions can be explained by examining the position- 
dependent LDOS 

N(z,E) = ~ J dk^MM^G^mhiz)), 

(4) 



from the left to the right Ni electrode. The LDOS is plot- 
ted in Fig. [5] where we choose four bias voltage values 
(V b = 0.0,-0.1,-0.5,-0.9 V) at which the magnitude 
of the total charge current differs significantly. In equi- 
librium (Vb = V), a prominent resonant state (white 
and red region) in the central Gr$ region is located close 
to the Fermi level and couples well with the conduction 
states on both sides of Ni electrode. Upon application 
of the bias voltage, both the position and the width of 
resonant states start to change. At Vb = —0.1 V, part of 
the resonant state enclosed by the dotted oval still follows 
rigidly the upward-moving conduction state in the left Ni 
electrode while extending all the way to the downward 
moving conduction state in the right Ni electrode. As a 
result of this strong coupling between resonant conduc- 
tion states within the energy interval [fiL^n] enclosed 
by the electrochemical potentials of the two electrodes, 
the current increases notably. However, the charge den- 
sity between the electrodes rapidly gets modified with 
the application of higher bias, and at Vb = —0.5 V, the 
resonant state splits into two parts (at energies ±0.2 eV) 
thereby losing coupling to one of the Ni electrodes. Thus, 
almost no current flows at this bias voltage. Increasing 
Vb further introduces a new state in the central region at 
energy w —0.4 eV which couples strongly to both Ni elec- 
trodes at Vb = — 0.9 V so that current starts increasing 
again. 



V. CONCLUDING REMARKS 

In conclusion, we demonstrated that perfect spin fil- 
tering in Ni|Gr n |Ni, with n > 5 layers of graphene sand- 
wiched between two (111) fee Ni electrodes, characte rized 
by "pessimistic" TMR=100% at zero bias voltage™^ 
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would persist even at finite bias voltage V& < 0.4 V. This 
feature is markedly different from conventional MgO- 
based MT Js where TMR drops sharpl^EEH with increas- 
ing bias voltage. Thus, it could play an important role 
in spintronic devices based on STTPEJ Furthermore, we 
predict that Ni|Gr n |Ni junction with P orientation of the 
Ni magnetizations would exhibit negative differential re- 
sistance as the bias voltage is increased from Vb = V 
to Vb — 0.5 V due to transmission resonance formed at 
zero bias voltage which is then gradually pushed outside 



of the bias window. 



ACKNOWLEDGMENTS 

We thank K. S. Novoselov and J. Q. Xiao for illumi- 
nating discussions. This work was supported by DOE 
Grant No. DE-FG02-07ER46374 (K. K. S. and B. K. 
N.) and Danish National Research Foundation's Center 
for Nanostructured Graphene (K. S. T.). The supercom- 
puting time was provided in part by the NSF through 
XSEDE resource TACC Ranger under Grant No. TG- 
DMR100002 and NSF Grant No. CNS-0958512. 



bnikolic@udel.edu 

1 W. H. Butler, X.-G. Zhang, T. C. Schulthess, and J. M. 
MacLaren, Phys. Rev. B 63, 054416 (2001). 

2 J. Katine and E. E. Fullerton, J. Magn. Magn. Mater. 320, 
1217 (2008). 

3 D. Ralph and M. Stiles, J. Magn. Magn. Mater. 320, 1190 

(2008) . 

4 C. Wang et al, Nature Phys. 7, 496 (2011). 

5 S. Yuasa, T. N. A. F. Y. Suzuki, and K. Ando, Nature 
Mater. 3, 868 (2004). 

6 J. Mathon and A. Umerski, Phys. Rev. B 74, 140404 

(2006) . 

7 Y. Ke, K. Xia, and H. Guo, Phys. Rev. Lett. 105, 236801 
(2010). 

8 P. Bose, A. Ernst, I. Mertig, and J. Henk, Phys. Rev. B 
78, 092403 (2008). 

9 D. Waldron, L. Liu, and H. Guo, Nanotech. 18, 424026 

(2007) . 

10 I. Rungger, O. Mryasov, and S. Sanvito, Phys. Rev. B 79, 
094414 (2009). 

11 V. M. Karpan et a/., Phys. Rev. Lett. 99, 176602 (2007). 

12 V. M. Karpan et al, Phys. Rev. B 78, 195419 (2008). 

13 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 

(2009) . 



14 C. Zhang et al, Phys. Rev. B 69, 134406 (2004). 

15 D. A. Areshkin and B. K. Nikolic, Phys. Rev. B 81, 155450 
(2010). 

16 T. Christen and M. Biittiker, Europhys. Lett. 35, 523 
(1996). 

17 J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 245407 
(2001). 

18 M. Brandbyge, J.-L. Mozos, P. Ordejon, J. Taylor, and K. 
Stokbro, Phys. Rev. B 65, 165401 (2002). 

19 Y.-H. Lee, Y.-J. Kim, and J.-H. Lee, Appl. Phys. Lett. 98, 
133112 (2011). 

20 M. Fuentes-Cabrera, M. I. Baskes, A. V. Melechko, and 
M. L. Simpson, Phys. Rev. B 77, 035405 (2008). 

21 M. Vanin et a/., Phys. Rev. B 81, 081408 (2010). 

22 F. Mittendorfer et al, Phys. Rev. B 84, 201401 (2011). 

23 T. Olsen, J. Yan, J. J. Mortensen, and K. S. Thygesen, 
Phys. Rev. Lett. 107, 156401 (2011). 

24 Y. S. Dedkov and M. Fonin, New J. Phys. 12, 125004 
(2010). 

25 J. Enkovaara et al, J. Phys.: Condens. Matter 22, 253202 
(2010). 

26 Atomistix ToolKit version 11.8, 
|http: / /www.quantumwise.com| 

27 M. A. Kuroda, J. TersofT, D. M. Newns, and G. J. Martyna, 
Nano Lett. 11, 3629 (2011). 



